UVIS Time-dependent photometry¶
Description¶
This jupyter notebook demonstrates how to apply the new time-dependent inverse sensitivities (zeropoints) using the PHOTFLAM header keyword to sample data.¶
The notebook comes with example standard star (GD153) flat-fielded CTE-corrected (FLC) files for three separate epochs in the filter F606W and C512C subarray mode. It also comes with a CSV file that contains a list of FLCs, and the centroids of the star for each FLC image.
Table of Contents¶
1. Import all modules¶
2. Read in CSV data, set up pandas dataframe, load pixel area map images¶
3. Compute Aperture photometry on the FLCs¶
4. Compute magnitudes, applying the new time-dependent PHOTFLAM (inverse sensitivity) keywords¶
5. Plot countrates and magnitudes¶
6. Equalize the countrate values in the science array of all FLC frames (optional step to prepare for drizzling).¶
7. Recompute aperture photometry and plot new countrates¶
8. Compute visit-level drizzled DRC data products¶
9. Recompute aperture photometry on the DRC files¶
10. Plot and compare all results¶
Imports¶
import glob
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from astropy.io import fits
from astropy.stats import sigma_clipped_stats
from astropy.coordinates import SkyCoord
from astropy import wcs
from drizzlepac import photeq
from drizzlepac import astrodrizzle
from photutils import aperture_photometry, CircularAperture, CircularAnnulus
from stwcs import updatewcs
%matplotlib inline
df = pd.read_csv('GD153_F606W_public.csv')
df
Pixel Area Maps (PAM)¶
FLC frames are not corrected for distortion and pixels therefore do not have equal area on the sky. To correct for this affect, we multiply the FLC frames by the Pixel Area Map. If you're using full-frame UVIS images, you do not need the dictionary below. Since the GD153 data we're looking at here are C512 subarrays, the PAM needs to be "cut out" at the region corresponding to the subarray.
pam1 = fits.getdata('UVIS1wfc3_map.fits')
pam2 = fits.getdata('UVIS2wfc3_map.fits')
pams = {'A': pam1[-512:, :513], 'B': pam1[-512:, -513:], 'C': pam2[:512, :513], 'D': pam2[:512, -513:]}
Photometry with no correction¶
Here we set up the code that computes the photometry. We are using a standard aperture size of 10 pixels, with a sky annulus from 155 to 165 pixels. These are the steps involved in computing the photometry:
- Loop through each FLC and load the data, MJD date of exposure, exposure time, and PHOTFLAM
- Divide the data by the exposure time to get pixel values in electrons per second (as opposed to just electrons)
- Multiply this result by the PAM, or the cutout of the PAM corresponding to the specific region of the subarray
- Supply the centroid of the object from the CSV file/pandas dataframe
- Set up the aperture, annulus aperture, and annulus mask
- Compute the sigma-clipped mean of the annulus
- Compute the photometry in the aperture
- Subtract the background measurement (derived in step 6) from the photometry derived in step 7.
- Store all values
phots = []
mjds = []
ap = 10
skyrad = [155, 165]
pfl = []
for i, flc in enumerate(df['FLC'].values):
with fits.open(flc) as f:
data = f[1].data
mjd = f[0].header['EXPSTART']
exptime = f[0].header['EXPTIME']
pfl.append(f[0].header['PHOTFLAM'])
data = data / exptime
data = data * pams[df.at[i, 'Amp']]
positions = (df.at[i, 'Centx'], df.at[i, 'Centy'])
aperture = CircularAperture(positions, ap)
annulus_aperture = CircularAnnulus(positions, r_in=skyrad[0], r_out=skyrad[1])
annulus_masks = annulus_aperture.to_mask(method='center')
annulus_data = annulus_masks.multiply(data)
mask = annulus_masks.data
annulus_data_1d = annulus_data[mask > 0]
mean_sigclip, _, _ = sigma_clipped_stats(annulus_data_1d)
background = mean_sigclip * aperture.area
apers = [aperture, annulus_aperture]
phot_table = aperture_photometry(data,apers)
final_sum = phot_table['aperture_sum_0'] - background
phots.append(final_sum[0])
mjds.append(mjd)
df['Countrate'] = phots
df['MJD'] = mjds
df['PHOTFLAM'] = pfl
Magnitudes¶
We can convert countrates into ST magnitudes using the zeropoint indicated in PHOTFLAM, and using the following equation. The EE_r10 is the encircled energy term for an aperture radius of r=10 pixels (0.4 arcseconds). For F606W, this is 0.91. This value can be computed using pysynphot.
EE_r10 = 0.91
df['STMags'] = -21.1 -2.5*np.log10(df['PHOTFLAM']) -2.5*np.log10(df['Countrate']) - 2.5 * np.log10(1./EE_r10)
df
Plotting Countrates¶
We first plot the photometric countrates (electrons per second) over time in MJD. Note that you can see the decline in the countrates over time, due to the expected and documented sensitivity losses for the instrument
fig = plt.figure(figsize=(40, 25), dpi=40)
plt.plot(df['MJD'], df['Countrate'], 'o', markersize=40, label='Amp C')
plt.grid()
plt.xlabel('MJD', fontsize=40)
plt.xticks(fontsize=40)
plt.yticks(fontsize=40)
plt.ylabel('Countrate (e-/s)', fontsize=40)
plt.title('GD153, F606W, UVIS2', fontsize=40)
plt.ylim(101000, 104000)
plt.legend(loc=0, fontsize=40)
Plotting Magnitudes¶
Now we plot the ST magnitude over time in MJD. This is using the corrected PHOTFLAM keyword, so as you can see the magnitudes are stable over time.
fig = plt.figure(figsize=(40, 25), dpi=40)
plt.plot(df['MJD'], df['STMags'], 'o', markersize=40, label='Amp C')
plt.grid()
plt.ticklabel_format(useOffset=False)
plt.xlabel('MJD', fontsize=40)
plt.xticks(fontsize=40)
plt.yticks(fontsize=40)
plt.ylabel('STMAG (Magnitude)', fontsize=40)
plt.ylim(13.55, 13.65)
plt.title('GD153, F606W, UVIS2', fontsize=40)
plt.legend(loc=0, fontsize=40)
Photometric equalization¶
This simple single-line step will equalize the countrates in the science array of the FLC frames to match any specified 'reference' image. For more details, see the phot_eq software documentation. Note that at this step we overwrite the science pixels in the original FLC files. Note that in this case the data are sorted, and the step automatically uses the 2009 GD153 data PHOTFLAM values to match the rest of the images to. You can supply a given reference PHOTFLAM value to the photeq call or/and ensure that your data are time-sorted
photeq.photeq(','.join(df['FLC'].values), readonly=False)
Repeating photometric Countrate calculation¶
We repeat the same computation as earlier for the photometric Countrate, but using the now photometrically equalized data.
phots = []
for i, flc in enumerate(df['FLC'].values):
with fits.open(flc) as f:
data = f[1].data
exptime = f[0].header['EXPTIME']
data = data / exptime
data = data * pams[df.at[i, 'Amp']]
positions = (df.at[i, 'Centx'], df.at[i, 'Centy'])
aperture = CircularAperture(positions, ap)
annulus_aperture = CircularAnnulus(positions, r_in=skyrad[0], r_out=skyrad[1])
annulus_masks = annulus_aperture.to_mask(method='center')
annulus_data = annulus_masks.multiply(data)
mask = annulus_masks.data
annulus_data_1d = annulus_data[mask > 0]
mean_sigclip, _, _ = sigma_clipped_stats(annulus_data_1d)
background = mean_sigclip * aperture.area
apers = [aperture, annulus_aperture]
phot_table = aperture_photometry(data,apers)
final_sum = phot_table['aperture_sum_0'] - background
phots.append(final_sum[0])
mjds.append(mjd)
df['Phot-eq'] = phots
df
Plotting corrected Countrates¶
We plot the photometric Countrate (electrons per second) over time in MJD. You can now see that the data are much less variant over time. This is the corrected data
fig = plt.figure(figsize=(40, 25), dpi=40)
plt.plot(df['MJD'], df['Countrate'], 's', markersize=50, label='Not equalized', alpha=0.55, color='Grey')
plt.plot(df['MJD'], df['Phot-eq'], 'o', markersize=40, label='Equalized')
plt.grid()
plt.xlabel('MJD', fontsize=40)
plt.xticks(fontsize=40)
plt.yticks(fontsize=40)
plt.ylabel('Countrate (e-/s)', fontsize=40)
plt.title('GD153, F606W, UVIS2 - Photometrically equalized', fontsize=40)
plt.ylim(101000, 104000)
plt.legend(loc=0, fontsize=40)
Performing Astrodrizzle and repeating photometry¶
The photometrically corrected data can be drizzled together (on a visit/epoch level basis) to correct for the effects of cosmic rays and to improve the signal-to-noise ratio.
We then perform the photometry on the finalized drizzle (DRC) products, using similar techniques as before. Drizzled data are in units of electrons per second and have the PAM already applied, so we no longer have to perform those two steps during photometry.
for i in range(max(df['Epoch'])):
flcs = list(df.query('Epoch == {}'.format(i+1))['FLC'].values)
astrodrizzle.AstroDrizzle(flcs,
output= 'allCframes_{}'.format(flcs[0].split('_flc.fits')[0]),
skymethod= 'match',
skystat='mean',
driz_sep_bits='80',
combine_type='median',
combine_nhigh=1,
driz_cr_snr= '3.5 3.0',
driz_cr_scale= '2.0 1.5',
final_bits= '80',
build=True,
clean=True,
preserve=False)
drcs = ['allCframes_ibcda4d0q_drc.fits', 'allCframes_ich304prq_drc.fits', 'allCframes_idbha6ntq_drc.fits']
drcents = [(327.9, 341.3), (251.2, 350), (266.6, 273.6)]
phots = []
mjds = []
pfls = []
for i, drc in enumerate(drcs):
data = fits.getdata(drc, ext=1)
mjds.append(np.mean(df.query('Epoch == {}'.format(i+1))['MJD']))
pfls.append(df.query('Epoch == {}'.format(i+1))['PHOTFLAM'].values[0])
positions = drcents[i]
aperture = CircularAperture(positions, ap)
annulus_aperture = CircularAnnulus(positions, r_in=skyrad[0], r_out=skyrad[1])
annulus_masks = annulus_aperture.to_mask(method='center')
annulus_data = annulus_masks.multiply(data)
mask = annulus_masks.data
annulus_data_1d = annulus_data[mask > 0]
mean_sigclip, _, _ = sigma_clipped_stats(annulus_data_1d)
background = mean_sigclip * aperture.area
apers = [aperture, annulus_aperture]
phot_table = aperture_photometry(data,apers)
final_sum = phot_table['aperture_sum_0'] - background
phots.append(final_sum[0])
mags = -21.1 - 2.5*np.log10(pfls[0]) - 2.5*np.log10(phots) - 2.5*np.log10(1./EE_r10)
Plotting¶
The corrected photometric Countrate from the DRC and the recomputed magnitudes are plotted below. Note how they agree with the results from the previous step using the FLC and the PAM.
fig = plt.figure(figsize=(40, 25), dpi=40)
plt.plot(df['MJD'], df['Phot-eq'], 'd', markersize=40, color='grey', label='FLC photometrically equalized', alpha=0.5)
plt.plot(mjds, phots, 'o', markersize=40, label='DRC photometry')
plt.grid()
plt.xlabel('MJD', fontsize=40)
plt.xticks(fontsize=40)
plt.yticks(fontsize=40)
plt.ylabel('Countrate (e-/s)', fontsize=40)
plt.title('GD153, F606W, UVIS2, DRC', fontsize=40)
plt.ylim(101000, 104000)
plt.legend(loc=4, fontsize=40)
fig = plt.figure(figsize=(40, 25), dpi=40)
plt.plot(df['MJD'], df['STMags'], 's', markersize=50, label='FLC', alpha=0.4, color='Grey')
plt.plot(mjds, mags, 'o', markersize=40, label='DRC')
plt.grid()
plt.ticklabel_format(useOffset=False)
plt.xlabel('MJD', fontsize=40)
plt.xticks(fontsize=40)
plt.yticks(fontsize=40)
plt.ylabel('STMAG (Magnitude)', fontsize=40)
plt.ylim(13.55, 13.65)
plt.title('GD153, F606W, UVIS2, DRC', fontsize=40)
plt.legend(loc=4, fontsize=40)
